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Abstract 

Four, linear, exponential, integration algorithms (two implicit, one explicit and one predic- 
tor/corrector) are applied to a viscoplastic model to assess their capabilities. Viscoplasticity 
comprises a system of coupled, nonlinear, stiff, first order, ordinary differential equations which 
are a challenge to integrate by any means. Two of the algorithms (the predictor/corrector and 
one of the implicits) give outstanding results, even for very large time steps. 
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is an increment in time. 

is the hardening modulus for yield strength. 

is the thermal diffusivity. 

is the bulk modulus. 

is the shear modulus. 

is the Kronecker delta; 1 if i — j, otherwise 0. 

is the infinitesimal strain tensor. 

is the (deviatoric) plastic strain tensor. 

is the Cauchy stress tensor. 

is the effective stress, 

is the magnitude of plastic strain rate. 

is the magnitude of effective stress. 

is the Macauley bracket of x\ x if x > 0, otherwise 0. 


2 Introduction 

In many scientific fields, systems of first order, ordinary differential equations are used to mod- 
el physical processes of interest. These equations are often coupled, mathematically stiff and/or 
nonlinear, thereby requiring special numerical algorithms for their solution [1]. Viscoplasticity is 
an example of such a modelling effort. The literature abounds with papers written on numeri- 
cal integration methods applicable to viscoplasticity [2, 3, 4, 5, 6, 7], all of which testify to its 
complexity. 

The paper begins with a presentation of four, linear, exponential, integration algorithms that 
have been recently derived by the authors [8]: two are implicit, the third is explicit, and the fourth is 
a predictor/corrector. The next section presents a viscoplastic model applicable to metals and solid 
solution alloys at high homologous temperatures. The final section presents numerical solutions for 
the viscoplastic model using the various integration algorithms. An appendix provides a discussion 
of Newton-Raphson iteration which we use in the solution schemes of the implicit integrators. 


3 First Order ODEs 

Viscoplasticity presents itself as a system of N (typically 3), coupled, nonlinear, stiff, first order, 
ordinary differential equations of the form 

Xa + U Q X Q = V Q (for a = 1 , 2 ,..., A), (1) 

where the X Q [t] are the N independent scalar, vector or tensor variables to be solved for. The 
dot 1 ‘ ’ is used to denote differentiation with respect to time, t. The square brackets [•] are used 
to denote ‘function of’, while parentheses (•) and curly brackets {•} are used for mathematical 
groupings. 

The parameters U Q [A^[t], tj and V 0 [X/ 3 [i],i are, in general, functions of the variables Xp and 
t. If neither parameter is a function of the independent variables, then the system of equations 
is said to be linear; otherwise, it is nonlinear — as is the case in viscoplasticity. To simplify the 
notation, we use square brackets containing time to indicate the dependence of any variable on 
both Xp[t] and t. For example, we write U Q [t + A/] and V a [t + At] to denote the values of the 

parameters U a J Xp[t + At],t + At and V a Xp[t + A t],t + At at time t + At. 
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3.1 Two Implicit Algorithms 

A technique like Newton-Raphson iteration (see the Appendix) must be used whenever a system of 
first order, ordinary differential equations is to be solved simultaneously using an implicit integration 
algorithm. We shall consider the following two implicit solutions for a system of N independent 
differential equations, viz . [8] 

X4t + A ,] = *„[«] + v a [, + A,] (- t , o|| + A , 1 , A1 ) A<, (2) 

for a = 1,2 , . . . , A, and 

X a [t + A t\ = X a [t] |C«-W+«-l*+ ac])-a* + 1 e -i(t/a[t]+f/o[t+At]).A( + y^ t + A< ]J . At , ( 3 ) 

for a = 1,2, . . N . These are the linear, asymptotic (2) and trapezoidal (3), implicit, exponential 
solutions for the integration of a system of ordinary differential equations (1). They are referred to 
as linear solutions becuase they were derived using Taylor and Euler-Maclaurin series expansions, 
respectively, that were truncated after their linear terms. Higher order (including exact) solutions 
can also be found in [8]. 

For large time steps, At >> 1, with U a > 0, the exponential term becomes small compared with 
1, and the asymptotic expansions 


lim X Q [t + At] x Vaf< + At \ 

At— ► large U a [t -f A<] 

0) 

lim X a [t + A<] x kV a [t + A<]-AZ 

Ai— ►large 1 

(5) 


are obtained for the asymptotic (2) and trapezoidal (3) solutions, respectively. The asymptotic 
expansion given in (4) is the correct asymptotic solution for the first order differential equation 
(1). Hence, there is an advantage in using (2) over (3) — the asymptotic solution (2) is accurate 
and stable for time steps of all sizes for exponentially decaying solutions, like those which occur in 
viscoplasticity. 


3.2 Explicit Algorithm 


Different needs demand different solution techniques, and Newton-Raphson iteration will not always 
be the technique of choice. Newton-Raphson iteration is ideally suited for large computer codes 
where large time steps are needed to minimize the computation time. But it is not the best choice 

when small time steps are required for gaining a detailed picture of a response. For applications of 

this type, explicit solution techniques usually work best. 

The linear, explicit, exponential, integration algorithm derived in [8] is given by 

+ Al) = + V Ui] (-■ £[,] A| I At, (6) 

for a = 1,2, . . ., N, Equation (6) differs from equation (2) in that all the terms on the right hand 

side of (6) are known explicitly; whereas, most of them in equation (2) are only known implicitly. 

For large time steps, At 1, with U Q > 0, the exponential term becomes small compared with 
1, and the asymptotic expansion 


lim 

At — Uarge 


A a [£ T Ai] 


Vo\A 

u a [t] 


(7) 
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is obtained for the explicit solution (6), Although the differences between the asymptotic (2) and 
explicit (6) solutions become negligible for small At, the differences between their asymptotic ex- 
pansions, equations (4) and (7), can be enormous for large At. If we start at t = 0 and apply a 
large time increment At, the solution V Q [At]/U a [At\ for the implicit approximation is asymptoti- 
cally correct, whilst that of the explicit solution, V a [Q]/U a [Q\, corresponds to a ratio of the slopes 
of V Q to U a at the initial time, t = 0. In the case of viscoplasticity, these situations correspond 
to a correct viscoplastic solution in the implicit approximation, and an incorrect elastic solution 
in the explicit approximation, respectively. In subsequent time steps the explicit approximation 
will oscillate around the true solution, but it will not become unstable. These oscillations can be 
mitigated only by choosing smaller time steps. 


3.3 Predictor/Corrector 


Predictor/corrector methods use an explicit algorithm to first look ahead and predict the value 
of the function at a future point. Knowing something about the future response of this function, 
one then uses an implicit algorithm to go back and make a second, more accurate, correction of 
what the future value of this function is. These methods have the advantage of being able to use 
the difference between predicted X' a [t + At] and corrected X"[t + At] values to establish an error 
which one can use to control the time step size, although that is not done herein. 

A predictor/corrector based on linear exponential solutions [8] can be constructed by considering 

/l_^aW-AA 

X' a [t + At] = X a [t] e-W A( + Va[t) f Ua[t] , At ) - At (for a = 1 , 2 ,..., TV) (8) 


for the predictor, and 

X^[t + At] = + V^[t + At] 


1 _ c -l/i[£+AtJ.At' 
Ujjt + A<]-A< 


•At (for a = 1,2, . ..,7V) (9) 


for the corrector, where 

U' a [t + At] = U a \X'p[t + At],t + At] and V^[t + At] = V a \x' fi [t + At],* + At 


( 10 ) 


The predictor (8) is underdamped while the corrector (9) is overdamped; therefore, an improved 
estimate for updating the variable X a is obtained by averaging or weighting their values, i.e . 

X a [t + At] = | + At] + X'a[t + A t]j . (11) 

In this integration algorithm, both the predictor and corrector are unconditionally stable, but the 
time step size must still be monitored for accuracy. 

When considering either the asymptotic, explicit or predictor/corrector algorithm in the neigh- 
borhood of {7 a -At « 0, one needs to expand (1 - exp[— U a 'At])/U a -At into a power series in order 
to secure a sound computational algorithm. 

These four, exponential, integration methods will now be applied to a viscoplastic model rep- 
resentative of copper at elevated temperatures. 


4 Viscoplastic Model 

Here we present a viscoplastic model whore the governing system of differential equations is suffi- 
ciently complex to warrant a Jacobian of dimension greater than 1. Our implicit integration algo- 
rithms have an important advantage over classical algorithms — like backward Euler integration — in 
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that they have the distinct possibility of substantially reducting the order of the Jacobian. In this 
viscoplastic model there is a reduction in order from a 13 X 13 to a 2 X 2 for the Jacobian matrix 
(13 is the spatial dimension for this system of equations). This reduction in order is of importance 
from a computational viewpoint. For a more detailed discussion of viscoplasticity in general, see 
Lemaitre and Chaboche [9]. 

Viscoplastic constitutive models are currently finding application in describing the rate depen- 
dent plastic response of metallic structures that undergo significant temperature change over their 
duty cycle. The isotropic, stress/ strain, constitutive equations are given by Hooke’s law, 


Sij = 2 n ( Eij - £?•) with e p kk = 0, 


( 12 ) 


°kk = 3k {e kk - a(T - T 0 )6 kk } , 


(13) 


where one sees that the plastic and thermal strains, £ P j and ot(T — 2o)^ij, are eigenstrains for the 
deviatoric and hydrostatic responses, respectively. The evolution of plastic strain is described by a 
flow law, 



(14) 


a kinetic law, 


||e p || = i?[T] Z 


r <M- y) i 


(15) 


and two evolutionary laws, 

Y = rj (VlP P ll - r[T,Y]), 


(16) 


(17) 

Viscoplasticity is an internal state variable theory where the internal state variables, Y and B tJ in 
this case, are governed by separate evolution equations. The evolution of internal state follows a 
competitive process between hardening and recovery mechanisms (both thermal and dynamic). The 
yield strength Y is a phenomenological representation of material strength; it reflects the density of 
dislocations. The back stress is a phenomenological representation of an internal stress state; 
it reflects the stress fields set up by dislocations in their heterogeneous substructures. 

We will consider a high-temperature viscoplastic model— valid for T > \T m , where T m is the 
absolute melting temperature — which is characterized by material functions that are power-laws; 
in particular: 


m = 

cxp [jf]’ 


\(\\Z\\-Y)} 


f<M- Y)\ 3 


D 


\ D ) ’ 


h[Y] = 

f Y \ 3 ~ n 
\yC) ’ 

(18) 
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Because of the chosen forms for the material functions h, L and r, this viscoplastic model analytically 
reduces to Odqvist’s [10] classical theory of creep at steady state, where there is no evolution of 
the internal state variables, and where 




defines the creep rate, with ||S|| = y\S{jSij characterizing the magnitude of deviatoric stress. 

For illustrative purposes, we consider the material constants given in Table 1, which approximate 
copper behavior in the neighborhood of 500° Cin the absence of dynamic recrystallization. 


5 Numerical Algorithms 


If we differentiate equation (12), combine it with (14), and subtract (17) from both sides; and then 
rearrange equations (16) and (17), the following set of differential equations is obtained: 


* . O' + W’iu _ „„p 

T n nil "t'j — + rrv -| 


Y +T) 


Z\\ ~ r ’~' 3 1 L[Y) 

r[r,r]-Mr]||e p lh v 


o, 


. H \ m \ R ^ ii* p ii ^ 

Uii + L[Y] Uij ~ \\E\\ ’ j ‘ 

This system of equations, N = 3, has the form of equation (1) where 

X x = Ey, X 2 = Y, X 3 = Bij , 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 


Ux = 






U 3 = 


i 

L[Y) ’ 


and 


* , m P \\ n 


V 2 = 0, 


V7 . r 


li^ll 


(24) 

(25) 


These are the governing differential equations for our viscoplastic model when cast into the form 
of equation (1). 

Upon examining equations (20-22), it becomes apparent why our Jacobian is at most a 3 x 3 
matrix (one for each of the three differential equations), whereas the Jacobian for backward Euler 
integration is a 13 X 13 matrix (one for each of the thirteen spatial dimensions: six for the effective 
stress, six for the back stress, and one for the yield strength). However, our Jacobian must be further 
reduced to a 2 X 2 matrix. This is because U\ and U 3 both become asymptotically proportional 
to ||e p || at steady state, and are therefore linearly dependent at steady state. Hence, either the 
evolution equation for back stress (22) or the evolution equation for effective stress (20) must not 
contribute to the construction of the Jacobian in order to prevent the Jacobian from becoming 
singular. Both approaches will work; however, they are not equivalent in computational efficiency. 
It is much more economical to construct the Jacobian matrix using equations (20 and 21) than 
using equations (21 and 22), because equation (20) contains information about the elastic response 
which is not present in equation (22). 
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5.1 Asymptotic Algorithm 

The linear, implicit, asymptotic, exponential, integration algorithm, equation (2), can be written 
as 

X a [t + A<] = X 0 [*]e~ ea + V a [t + AZ] ( — ^ — j At for a = 1,2 (26) 

where, for our viscoplastic model, 


q. + ff)P-[ 
Sl IIEII 


and 


Q2 


= 7? ( 


r\TX\- h[Y}\\e”\ 


O’ 


which are both evaluated at time t + At. Thus, the iteration functions become 


(27) 




\\z[Mx}\\ 




(28) 


and 


^2 [{gp} a ] 


V (r| 


- h 

Y[{g 2 } x] 

P P [to,f?2}A]||) 

Y 

[{S2}X) 


{Q2}\, 


(29) 


which we have solved through Newton- Raphson iteration (c/. Appendix). The derivatives required 
to construct this 2x2 Jacobian matrix, e.g. d\\e p \\ jdgx, etc have been determined numerically, 
for in this case the determination of numerical derivatives is computationally more efficient than 
determining them analytically. 

The Newton-Raphson iterations are accomplished by the following scheme. First, values of Q\ 
and g 2 are guessed. We set their values to 0.1 for the first time step, and used the converged values 
of the prior time step as the initial guess for all time steps thereafter. In addition, we also guess 
the values of Y and B t j in like manner. The forcing vectors, V 2 and V 3 are then computed 
and used with £1 and g 2 in the recursion relationships to determine Xi, X 2 and X 3 . Estimates 
of E,y, Y and Bij are now available to compute improved estimates of the forcing vectors, V\, V 2 
and V 3 for the next iteration. It is evident from the preceding that this algorithm is basically one 
iteration in arrears in estimating £ ZJ , Y and B XJ . 

The capability of this integration method is demonstrated in Fig. 1. The solid curve was 
obtained using 500 integration steps, and is considered to be the converged solution. In this 
example, a homogeneous block of material is sheared in a specified direction at a constant rate of 
straining to a fixed value of engineering strain, i.e . 7 = 2£ U . This block of material is then sheared 
in the opposite direction to a self-similar value of engineering strain. The direction of loading is 
changed one last time to complete the loading cycle. The differences between the converged response 
and those determined using 25 and 3 integration steps are very small; however, there is a small 
amount of measurable error for the case of 10 integration steps at the knee of the curve. 1 These 
errors are enumerated in Table 2. At first glance, the fact that 3 integration steps do better than 
10 seems contradictory. However the contradiction is apparent only, because the three integration 
points are close to their asymptotic solutions where the implicit, asymptotic, integration method 
is very accurate. The regions within which the 10 integration points are in greatest error are the 
transient domains where the back stress is rapidly evolving. 

One important observation about this integration method is that the error generated in the 
transient domains does not propagate with the solution into the asymptotic domains. This is 
because the correct asymptotic solution (4) of the differential equation (1) is contained within the 
linear, implicit, asymptotic, integration method. 


^or all of the results presented in this paper, none of the time 


steps were subincremented. 
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The increase in the shear stress from the first reversal to the final state is due to a gradual 
increase in the value of the yield strength Y over the loading history, i.e. the material is gradually 
getting stronger because of work hardening. The smooth curvature observed in the stress-strain 
response of each of the three loading segments is due to the rapid evolution of the back stress Bij 
over each segment. 

5.2 Trapezoidal Algorithm 

For implicit, trapezoidal, Euler-Maclaurin integration, one can write the exponential solution (3) 
as 


X a [t + At] = X a [t] e~ ea At + \ ( V a [t]e~ ea - At + V a [t + At]) -At for a = 1, 2 
where, for our viscoplastic model, 

L ±H_ ( \\i’[t]\\ ||g p [< + A<]|| \ 

“ 2 l||S[t]|| + \\Z[t + At}\\)' 

and 

2 ( r W ~ ^MII^MH , + Afl - + At]||e p [< + Ai]|| \ 

2 \ Y[t] + Y[t + At] ) ‘ 

The associated iteration functions are therefore given by 


[{^}a] = 


fi+H f\\^[t}\\ ||e p [{p 1 ,p 2 } A ]||\ 

2 VII^WH + ||^[{p 1 } A ]|| ) 


and 


(30) 


(31) 


(32) 


(33) 


•^2 [{ 0 /?} a ] 


l( r[t]-h[t] ||6 p [t]|| 
2 y Y[t } 


+ 


r.y[{g2}A]]-fe[y[{g 2 } A ]] iie p Kei, *>2}a] 


Y[{Q2}x] 


{Q2}\, 


(34) 


which are solved through Newton-Raphson iteration, as described in the Appendix. 

The capability of this integration method is shown in Fig. 2 for the same loading history that 
was used in Fig. 1. The solid curve represents 500 integration points, and is equivalent to the 
converged response of Fig. 1, except at the locations of load reversal. The errors incurred by using 
fewer time steps, viz, 25, 10 and 3, are given in Table 2. Clearly one must monitor the size of 
the time step when using trapezoidal integration so as to secure accurate answers. Of the four, 
exponential, integration methods investigated herein for the numerical integration of viscoplastic 
models, the trapezoidal method is the least desirable. However for other applications, it can be the 
preferred method of integration [8], 

It is apparent in Fig. 2 that the trapezoidal solution does not converge towards the asymptotic 
solution, as is the case with the asymptotic integration method. The inaccuracies attendant on 
the larger time steps are due to the representation in equation (20). The right hand side of this 
equation contains the term ||e p || which depends in a highly nonlinear manner on the variable 
on the left hand side of the equation. 
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5.3 Explicit Algorithm 


The linear, explicit, exponential, integration algorithm is much simpler than the prior two algo- 
rithms in that an iterative solution procedure is not required, and therefore there is no need to 
construct a Jacobian. The integration is effected by directly substituting the X’s, tTs and V’s of 
equations (23-25) into the recursive solution (6). 

The capability of this integration method is demonstrated in Fig. 3. Again, the solid curve 
represents the converged response. There is less error with this method than with the trapezoidal 
method for the larger time steps, as shown in Table 2. This is suprising at first glance because the 
trapezoidal method is implicit. The reason why the explicit algorithm does so well is that it contains 
the correct form for the asymptotic solution — the trapezoidal algorithm does not. Although the 
form of the asymptotic solution is correct, i.e. V a /U ai it is evaluated at time t instead of time 
t + A/, as it should be. As a consequence, several time steps must be incurred along the loading 
path to secure an accurate result, as observed in Fig. 3. If the time steps become too large the 
solution will oscillate, and the extent of oscillation will increase with the size of the time step. 

For comparative purposes with Fig. 3, calculations obtained with the well-known forward Euler 
method are presented in Fig. 4. Our explicit method is the exponential analog of the forward Euler 
method. A comparison of errors between the explicit and forward Euler methods is also given in 
Table 2. It is apparent from both the figures and the table that the explicit exponential method is 
better than the explicit Euler method. 


5.4 Predictor/ Corrector 

Like the explicit algorithm, the exponential predictor/corrector is simple to construct because a 
Jacobian is not required. The integration of the predictor follows the same logic as the explicit 
algorithm which was just discussed. With the predicted quantities X' a [t + At] now known, the (/’ s 
and V's are updated to time t + A t, and integration via the corrector is effected by substituting 
equations (23-25) into the recursive implicit solution (9). For the results presented herein, the 
predictor and corrector were both integrated once at each time step — there was no iteration of the 
corrector. 

The ability of this integration method is demonstrated in Figs. 5 and 6. As before, the solid 
curves represent the converged solution and were obtained using 500 integration steps. In Fig. 5, 
the results that are presented are the corrected (not weighted) values X ff coming from equation (9). 
In Fig. 6, the results that are presented are the averaged (weighted) values -(A' + X n ) coming from 
equation (11). In both cases, the presence of the corrector dampens the oscillations that are present 
in the explicit (predictor) solution, cf Fig. 3. The reason why an improved result is obtained by 
averaging the predicted and corrected values is because the predictor is underdamped while the 
corrector is overdamped. Consequently, averaging them tends to remove much of the unwanted 
damping. 

For comparative purposes with our weighted predictor/corrector, i.e . Fig. 6, calculations ob- 
tained with the well-known Heun method are presented in Fig. 7. A comparison of errors between 
these two methods is also given in Table 2. The weighted predictor/corrector is the exponential 
analog of Heun’s method, which is a second order Runge-Kutta method. It is apparent from both 
the figures and the table that the exponential, weighted, predictor/corrector is better than the 
Heun method. 


9 



6 Concluding Remarks 

Four, linear, exponential, integration algorithms that we derived in [8] have been applied to a 
viscoplastic model which is composed of a system of coupled, nonlinear, stiff, first order, ordinary 
differential equations. Of these four methods, two have been found to have superior properties for 
the solution of these equations. The implicit asymptotic algorithm is the most accurate, but it 
requires the construction of a Jacobian matrix. It is recommended for large computer codes like 
finite elements. The averaged or weighted predictor/corrector also has exceptional accuracy, plus 
it has the advantage of not needing a Jacobian. It is recommended for smaller codes, and for those 
cases when implicit algorithms are not practical, e.g. our composite micromechanics theory [11]. 

For exponentially decaying solutions, like those which occur in viscoplasticity, these two inte- 
gration algorithms have the following desirable properties: 

• They are asymptotically correct. 

• They are stable and accurate for both small and large time steps. 

• They do not propagate error. 

• They are computationally efficient. 

When compared with equivalent classical methods, the exponential integration algorithms are supe- 
rior. This is because viscoplasticity has exponential solutions for which the exponential integration 
algorithms are better suited. 
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Appendix: Newton-Raphson Iteration 


In general the g a , where g a = U a [t + At] in (2) or g a = ^(f/ a [<] + U a [t + At]) in (3), are 
functions of the X@[t + At] which are in turn, via (2) or (3), functions of the Q„\t + At], and 
therefore 


Qa 


X0 [ft,]] = Qa [f?i/] 


(for a,/3, v — 1 to AT), 


(Al) 


or equivalently, 


where 


F a [eA = o, 

F a [Qu\ — Qa [Qv\ Q<y 


(A2) 

(A3) 


If {^a}^ ls th.e A th guess for vector g Q , then the true solution satisfying (2) or (3) may be written 
as 


Qa — F 


(A4) 


where the correction vector, c a , is the amount by which the true vector differs from the guessed 
value. Inserting this definition into (A2), and expanding the resultant by Taylor’s theorem while 
retaining only the first order terms leads to 


= *.[{<>,}» + c,] a + E Cg = 0 . 

i3- 1 

A solution to the resulting linear system of equations, 


N 

^ Ja(3 c {3 — ~~F o [ { } A ] 
0=1 


(for a = 1 to N), 


(A5) 


(A6) 


for the unknown correction vector, cp, with Jacobian, 


_ gW a] 
HQffh ’ 


(A7) 


can be obtained by Cramer’s rule ( N < 3) or Gaussian elimination ( N > 4). This leads to an 
improved value for the solution vector, {p a },\+i> where 


{PoJa+I “ {Q«h + C on 


(A8) 
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which is iterated until the contribution from c a becomes negligible (we used ||c|| < 0.0001||^||a as 
our convergence criterion). 

In order to retain algorithmic stability, it is necessary to bound the size of the correction vector 


so that 

IMI < T Ilf ||a 

(A9) 

(we typically set T % 0.1) where 



||c|| = \/ C l 2 + C2 2 + --- + Cyv 2 

and Weh = \/ Wa + Mi + • • • + {9n)1 ■ 

(A10) 

That is, if 

set c a => c q T J# V a G {1, N] 

ll c ll 


Ikll > T II 0 II 

(All) 


in equation (A8). This slows down the rate of convergence, but it keeps the procedure stable 
and free of oscillations. What we are actually doing is keeping the solution within the domain of 
applicability of the truncated Taylor expansion given in equation (A5). 

If possible, the Jacobian (A7) should be determined analytically; if not, it can be acquired 
numerically, but often times at a greater numerical expense. Because Newton-Raphson iteration 
has quadratic convergence, the number of iterations A required for convergence is usually few in 
number. One reasonable initial guess to advance the iteration process to the next step is to use the 
values from the last time step, z.e. set g a [t -f A/]a=i = Qot[t]- 
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Table 1: Elastic and viscoplastic 


material constants. 


Constant 

Units 

Value 

a 

1/°C 

20 X 10~ e 

K 

MPa 

95,000 


MPa 

30,000 

c 


0.8 

D 

I 

0.016 

H 

l 

15,000 

n 


5 

Q 

J/mol. 

200,000 

y 


0.1 

V 


30,000 


Table 2: Percent error in the range of stress, 100-|(<7 - asooV^maxI, evaluated at two points in the 
third branch of each curve. Errors greater than 100% are denoted as >100. There the 
algorithm is outside its domain of applicability and the error is meaningless. 


Algorithm 

Type 

Integration Points 

£ = 

0.0 

£ 

= 0.005 


10 

25 

3 

10 

25 

Asymptotic 

9.8 

2.1 

1.8 

1.2 

0.3 

Trapezoidal 

26.2 

2.3 

>100 

38.7 

6.2 

Explicit 

19.6 

0.9 

56.5 

12.5 

0.0 

Forward Euler 

17.3 

19.8 

>100 

>100 

1.5 

Predictor/Corrector 

23.9 

5.0 

>100 

16.1 

0.7 

Weighted Pred/Corr 

1.3 

0.0 

45.3 

0.6 

0.3 

Heun 

62.0 

15.8 

>100 

26.0 

1.6 
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Shear Stress, MPa Shear Stress, MPa 



Figure 1. — Linear, implicit, asymptotic solu- 
tions for the viscoplastic representation of 
copper. 7 = ±0.001s -1 . T = 500 °C. Y| t = 0 
= 1 MPa. The integration points were evenly 
spaced: — represents 500 steps, x repre- 
sents 25 steps, + represents 10 steps, and 
o represents 3 steps. 



Figure 3. — Linear explicit solutions for the visco- 
plastic representation of copper. 7 = ± 0.001 s" 1 . 
T = 500 °C. Y| t = o = 1 MPa. The integration 
points were evenly spaced: — represents 500 
steps, x represents 25 steps, and + repre- 
sents 1 0 steps. Stresses for 3 integration steps 
exceeded the plotted range of stress. 



Figure 2. — Trapezoidal Euler-Maclaurin solu- 
tions for the viscoplastic representation of 
copper. 7 = ±0.001 s - ^ . T = 500 °C. Y| t = 0 
= 1 MPa. The integration points were evenly 
spaced: — represents 500 steps, x repre- 
sents 25 steps, and + represents 10 steps. 
Stresses for 3 integration steps exceeded the 
plotted range of stress. 



Engineering Shear Strain 

Figure 4 —Forward Euler solutions for the visco- 
plastic representation of copper. 7 = ± 0.001 s _1 
T = 500 °C. Y[ t = Q = 1 MPa. The integration 
points were evenly spaced: — represents 500 
steps and x represents 25 steps. Stresses for 
10 and 3 integration steps exceeded the plotted 
range of stress. 
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Engineering Shear Strain 

Figure 5. — Linear predictor/corrector solutions 
(not weighted) tor the viscoplastic represen- 
tation of copper. y - ±0.001 s -1 . T = 500°C. 
Y| t = o = 1 MPa. The integration points were 
evenly spaced. — represents 500 steps, x 
represents 25 steps, and + represents 10 
steps. Stresses for 3 integration steps ex- 
ceeded the plotted range of stress. 



Engineering Shear Strain 

Figure 6. — Linear predictor/corrector solutions 
(weighted) for the viscoplastic representation 
of copper. 7 = ±0.001 s -1 . T =500 °C. Y| t = 0 
= 1 MPa. The integration points were evenly 
spaced. — represents 500 steps, x repre- 
sents 25 steps, and + represents 10 steps. 
Stresses for 3 integration steps exceeded the 
plotted range of stress. 



Engineering Shear Strain 

Figure 7. — Heun (a second order Runge-Kutta) 
solutions for the viscoplastic representation of 
copper. 7 = ±0.00fs“ 1 . T = 500 °C. Y| t = 0 
= 1 MPa. The integration points were evenly 
spaced. — represents 500 steps and x re- 
presents 25 steps. Stresses for the 10 and 3 
integration steps exceeded the plotted range 
of stress. 
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